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An important part of building mathematical models based on measured data is calculating the accuracy associ- 
ated with statistical estimates of the model parameters. Indeed, without some idea of this accuracy, the parameter 
estimates themselves have limited value. An expression is developed for computing quantitatively correct param- 
eter accuracy measures for maximum likelihood parameter estimates when the output residuals are colored. This 
result is important because experience in analyzing flight test data reveals that the output residuals from maximum 
likelihood estimation are almost always colored. The calculations involved can be appended to conventional maxi- 
mum likelihood estimation algorithms. Monte Carlo simulation runs were used to show that parameter accuracy 
measures from the new technique accurately reflect the quality of the parameter estimates from maximum likeli- 
hood estimation without the need for correction factors or frequency domain analysis of the output residuals. The 
technique was applied to flight test data from repeated maneuvers flown on the F-18 High Alpha Research Vehicle. 
As in the simulated cases, parameter accuracy measures from the new technique were in agreement with the scatter 
in the parameter estimates from repeated maneuvers, whereas conventional parameter accuracy measures were 
optimistic. 


Nomenclature 

a z - vertical acceleration, g 

Cl = lift coefficient 

Cm = pitching moment coefficient 

Cz = vertical force coefficient 

c = mean aerodynamic chord, ft 

D = dispersion matrix 

djj = jth diagonal element of D 

E{ ] = expected value 

e(i) = ith equation error residual 

g - gravitational acceleration, 32.174 ft/s 2 

I y = pitch axis moment of inertia, slug-fit 2 

J - cost function 

M = information matrix 

m = mass, slug 

N = total number of sample times 

n 0 - number of outputs 

n p = number of parameters 

q = body axis pitch rate, rad/s 

q - dynamic pressure, lbf/ft 2 

R = discrete noise covariance matrix 

9^ = autocorrelation matrix of vector v 

S - wing area, ft 2 

S(i) = output sensitivity matrix at time (i — l)Af 
s = sample standard error 

t = time, s 

k( 0 = control vector 

V = airspeed, ft/s 

v(i) = output residual vector at time (i — l)Ar 
x(f) = state vector 

y(i) = n 0 xl output vector at time (i — 1) At 
y(t) = n 0 x 1 output vector 
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z(i) 

= ith measured scalar 

z(0 

= measured n a x 1 output vector at time (i - 1) At 

a 

= angle of attack, rad 

At 

= sample time, s 

&ij 

= Kronecker delta 

Ss 

= stabilator deflection, rad 

© 

= pitch angle, rad 

0 

= element of the parameter vector 6 

0 

= n p x 1 parameter vector 

<j 

= Cram£r-Rao bound for the standard error of 0 

v(i) 

= measurement noise vector at time (i - 1) At 

v* 

= gradient with respect to 8 

<l> 

= roll angle, rad 

0 

= zero vector 

Subscripts 

c 

= corrected 

m 

= measured 

o 

= initial or bias 

w 

= wind axes 

Superscripts 

T 

= transpose 

-1 

= matrix inverse 

* 

= estimate 


= mean value 


= time derivative 


Introduction 

A IRCRAFT dynamic models include parameters that quantify 
the dependence of aerodynamic forces and moments on state 
and control variables. The values of these parameters are often es- 
timated from flight test data. A good quantitative assessment of the 
accuracy of these parameter estimates is important for a variety 
of reasons related to experiment design, modeling, simulation, and 
flight control. 

Maximum likelihood 1,2 is commonly used to estimate aero- 
dynamic parameters from flight test data. Assuming the model 
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structure is correct, maximum likelihood parameter estimates ap- 
proach the true parameter values, and the parameter variances ap- 
proach their theoretical minimum values (the Cramfr-Rao lower 
bounds), as the number of measured data points increases. Gener- 
ally, a flight test data record length at least two to three times the 
period of the slowest dynamic mode to be modeled is sufficient 
for the parameter variances to closely approach the Cram6r-Rao 
bounds. 3 In such cases, the Cram6r-Rao bound can be used as a 
good approximation to the variance of maximum likelihood parame- 
ter estimates. References 3-5 compare and contrast the Cram£r-Rao 
bound with other methods for assessing the accuracy of parameter 
estimates. Theoretical properties of maximum likelihood estimators 
and related arguments discussed in Ref. 3 indicate that the Cram^r- 
Rao bound is the best accuracy measure for maximum likelihood 
parameter estimates. 

The research described here focuses on the output error formula- 
tion of maximum likelihood parameter estimation. This formulation 
includes measurement noise, but no process noise. 1 - 2 A modified 
Newton-Raphson optimization procedure 6 was used to determine 
the maximum likelihood parameter estimates. With this approach, 
the Cram6r-Rao bounds are computed as part of the estimation pro- 
cedure. It is well known, however, that the Cram6r-Rao bounds 
computed in this way are usually optimistic (too small) compared 
to the scatter in the parameter estimates from repeated flight test 
maneuvers. 3,7 This prompted the work of Maine and Diff 2 3 (also 
Refs. 8 and 9) and Balakrishnan and Maine, 10 who traced the dis- 
crepancy to the fact that the output residuals are colored for real 
flight test data analysis because some deterministic modeling error 
is always present. Output error techniques lump the deterministic 
modeling error together with the broadband random part of a mea- 
sured signal and call this the measurement noise. This means the 
measurement noise is model dependent and colored, because the de- 
terministic modeling error usually lies in the same frequency band as 
the aircraft rigid body dynamics and accounts for a large part of the 
total noise power. References 2, 3, and 8-10 describe how this kind 
of colored measurement noise is responsible for the discrepancy be- 
tween the conventional calculation of the Cram^r-Rao bounds and 
the observed scatter in flight-determined parameter estimates from 
repeated maneuvers. 

The theory underlying the output error formulation of maximum 
likelihood estimation assumes that the measurement noise is white 
Gaussian and band limited by the Nyquist frequency. The band 
limit is the result of discrete measurements taken at the sampling 
frequency, which is twice the Nyquist frequency. This measurement 
noise is broadband and incoherent. The term incoherent implies 
amplitude discontinuity and a lack of consistent phase-amplitude 
relationships, causing the autocorrelation function to be close to 
the impulse function. This part of the residual would be commonly 
recognized as having no deterministic component If the structure 
of the model were correct, the residuals would be expected to be 
reasonably close to this type of noise. In real flight test data analy- 
sis, however, the residuals contain deterministic components from 
such sources as approximations to real aircraft aerodynamic de- 
pendencies, unmodeled dynamics such as structural modes, and lin- 
earization of the equations of motion. The result is colored residuals, 
which violate the assumption of white measurement noise in conven- 
tional maximum likelihood theory, leading to the aforementioned 
discrepancy. 2,3, 8_ 10 

In Ref. 3, several engineering solutions were proposed to correct 
for the discrepancy. Each solution was based on the assumption that 
most of the power in the output residuals for real flight data analysis 
is concentrated in roughly the same frequency band as the rigid body 
dynamics and is due to deterministic modeling error. This assump- 
tion is stretched when relatively high-frequency structural modes 
appear in the data or when the broadband random noise has a large 
enough magnitude to rival the power of the narrow-band noise due to 
deterministic modeling error. For multiple outputs, the noise power 
from broadband random noise compared to that from narrow-band 
deterministic modeling error is different for each output because of 
differences in the sensor characteristics and the physical quantity 
being measured. The solutions offered in Ref. 3 depend on knowing 
something about the bandwidth of the dominant source of power 


in the residuals. Obtaining this information requires Fourier trans- 
forms of the residuals and analysis in the frequency domain. The 
spectra of the residuals depend on the model structure, the maneu- 
ver, the flight condition, and the instrumentation characteristics. All 
of these factors can change over the course of a flight test program, 
requiring changes in the corrections for the Cram6r-Rao bounds. 
In addition, the band widths of the deterministic modeling error for 
the various measured outputs can be different from one another for 
the same maneuver. The solutions from Ref. 3 require some engi- 
neering judgment in the form of determining a correction factor or 
estimating the bandwidth of the dominant power in the residuals. 
Both of these approaches require an experienced analyst and limit 
the accuracy of the results. 

In the present work, a technique first put forth in Ref. 1 1 was 
used to process the residuals from a conventional maximum like- 
lihood estimation to compute accurate Cram6r-Rao lower bounds 
for colored residuals. The approach accounts for colored residu- 
als using a simple estimate of the residual correlation in the time 
domain. Existing maximum likelihood estimation routines can be 
easily upgraded because the technique involves a postprocessing 
of the output residuals to correct the Cram6r-Rao lower bounds 
from the conventional calculation. The purpose of the present work 
is to document a few refinements and extensions, to validate the 
technique using Monte Carlo simulation with colored measurement 
noise, and to apply the technique to real data from repeated flight test 
maneuvers. 

The next section contains the theoretical analysis. Following this, 
the technique was applied in a controlled situation using simu- 
lated data from a model of the longitudinal dynamics of a fighter 
aircraft. The true parameter values were known, and the mea- 
sured outputs were corrupted with colored noise, including both 
narrow-band modeling error and broadband random noise. Using 
200 Monte Carlo simulation runs with various colored noise char- 
acteristics, it was demonstrated that the new technique produces 
Cram6r-Rao bounds representative of the observed scatter in the 
parameter estimates. The conventional Cramdr-Rao bounds were 
found to be optimistic, in agreement with the results of previous 
research. 2,3,7 " 10 

The technique was then applied to repeated longitudinal flight 
test maneuvers at 20-deg angle of attack for the F- 18 High Alpha 
Research Vehicle (HARV). The scatter in the model parameter esti- 
mates from this flight test data was consistent with the Cram6r-Rao 
bounds computed using the new technique, whereas the conven- 
tional calculation again gave optimistic values for the Cram6r-Rao 
bounds. 

Theoretical Development 

The aircraft dynamic model can be represented as 

*(0 = /W0.«(0.fl (1) 

x(0)=x o (2) 

y(t) = g[x(t\u(tie] 0) 

Z(i) = y(i) + v(i) i = 1,2 N (4) 

For conventional maximum likelihood, the discrete measurement 
noise vector u(i) is assumed to be zero mean white Gaussian and 
band limited at the Nyquist frequency, 

£{«(«)} = 0 E{v(i)v T (J)} = RS,j (5) 

The maximum likelihood estimate of the parameter vector max- 
imizes the conditional probability density function 1,6 

6 = argmax[p(Z | 0)] (6) 

where Z is the set of all measurement vectors z(i), for i — 
1,2, ...,V. When R is known, maximizing the conditional 
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probability in Eq. (6) is equivalent to minimizing the cost func- 
tion 1 ' 3 ' 6 ' 11 

1 N 

m = - £ wo -y(of* -1 [z(o ->(<)] (7) 


For practical computation, simultaneous satisfaction of the fol- 
lowing numerical criteria was used to define convergence of the 
maximum likelihood estimation: 

I [0y]* - [$,]*- 1 1 < 1.0 X 10- 5 V/\ ; = 1, 2 n„ 


The cost in Eq. (7) can be minimized using a modified Newton- 
Raphson technique 6 to determine parameter updates, starting from 
some initial guess of the parameter vector. The initial guess for the 
parameter vector can be obtained from equation error methods, 7 but, 
typically, a much rougher initial guess can be used. 

The sensitivity matrix is defined as 


Inih-ln r]*-i 
Ril*. i 


<0.05 


Vi, i = 1, 2, . . . , n 0 

( 12 ) 


j{e k )-j(e k . x ) 

J( 0 t- 1) 


< 0.001 


S(i)E 


By( 0 


d0 


10 = 0 


i = l,2 N (8) 


dj(9) 


dOi 


<0.05 


$=e 


j — 1, 2, . . . , Tip 


where the jth column of the sensitivity matrix contains the output 
sensitivities for the jth parameter, computed from central finite dif- 
ferences in Eqs. (1-3). The modified Newton-Raphson parameter 
update is given by 1,6,11 

ao = e-e 



S(f) r jr‘s(o 


N 

^SdfR-'izd) -yd)) 

/«! 


(9) 


The parameter vector update from Eq. (9) is added to the current 
estimate of the parameter vector to approach the true value of the 
parameter vector. In practice, there are times when the parameter 
vector update computed from Eq. (9) leads to an increase in the cost 
function or a divergence. This is because the modified Newton- 
Raphson step assumes that the current estimate of the parameter 
vector is near the true value. Using several iterations of a simplex 
algorithm 12 when the modified Newton-Raphson step produced an 
increase in the cost was found to be very effective in avoiding diver- 
gence and reaching a solution. This approach was followed in the 
present study. 

When repeated application of Eq. (9) converges, an estimate of 
the measurement noise covariance matrix R can be obtained from 
the output residuals. The expression for the estimate of R , which 
maximizes the conditional probability in Eq. (6), is 1-3,6 * 11 


R = ^ J^lzd) -ymizd) -yd)) T do) 


The most recent estimated output y(i ) is substituted for y(i) in 
Eq. (10) to compute the matrix R . Often only the diagonal elements 
of the R matrix are estimated from Eq. (10), enforcing an assump- 
tion that the measurement noise sequences for the measured outputs 
are uncorrelated with one another. This assumption is generally a 
good one for real flight test data. All estimates of the measurement 
noise covariance matrix in this work assume a diagonal R matrix. 
Retaining the full R matrix could have been done with little con- 
ceptual difficulty, but the expected benefits did not warrant the extra 
computation involved. 

The noise covariance matrix estimate R was used in the cost func- 
tion of Eq. (7), and the minimization process described earlier for 
known R was repeated. Thus, the maximum likelihood estimation 
proceeds by alternately estimating the noise covariance matrix from 
Eq. (10), and minimizing the cost function using Eq. (9) with the 
latest value of the estimated noise covariance matrix. Convergence 
is reached when the estimated parameter vector 0 , the estimated 
noise covariance matrix R, and the cost J(0 ) reach nearly constant 
values. 

Since maximum likelihood estimation^ is asymptotically un- 
biased, 1 ' 3 the estimated parameter vector 6 , should be close to the 
true value 6 , and the gradient of the cost function with respect to the 
parameter vector should be close to zero. From Eq. (7), assuming R 
is held fixed. 


V0/(0)l* = * 


-^S(«) r *-‘[z(0 ->(«)] 


i = I 


where k denotes the current estimate iteration number and ru de- 
notes the estimate of the ith diagonal element of R. The approxi- 
mate expression for the cost gradient with respect to the parameters 
[Eq. (1 1)] was used for the last criterion in Eq. (12). 

The minimum achievable parameter variances, called the Cram- 
6r-Rao lower bounds, are the diagonal elements of the dispersion 
matrix D (Refs. 1-3 and 6). The dispersion matrix is defined as 
the inverse of the information matrix Af, the latter being a measure 
of the information contained in the data from an experiment. The 
expressions for these matrices are 1 ' 3,6 

N 

M=^Sd) T R' l Sd) (13) 

i = 1 

d = m~ 1 = |^S(,yjr*S(0 

The square root of the jth diagonal element of D gives the 
Cramer-Rao lower bound for the standard error of the j th parameter 
estimate, 

= y = l,2, (15) 

It can be seen from Eqs. (9) and (14) that the dispersion ma- 
trix is computed when determining the modified Newton-Raphson 
step as part of the conventional maximum likelihood estimation. 
The assumption that the output residuals are white and, therefore, 
unconelated in time is implicit in the algorithm and indicated in 
Eq. (5). The next section details the theory involved in accounting 
for arbitrary colored output residuals, which are correlated in time. 

When the maximum likelihood estimation has converged, the 
estimated parameter vector will be close to the true value and Eq. (9) 
holds. Define the residual vector 

v(i) £z(i) -y(0 i = l,2, ...,JV (16) 

The estimated parameter covariance matrix can be expressed using 
Eq. (9) with substitutions from the definitions in Eqs. (14) and (16), 

cov(0) == E[(0 — 0){0 - Of] 

= e| Y^DSd) T R- l vd)v(J) T R~ l S(J)D j (17) 


(14) 


If it is assumed that the discrete noise covariance matrix and the 
output sensitivities do not depend on the parameter vector estimate 
at the maximum likelihood solution, then the estimated parameter 
covariance matrix can be written as 


cov(0 ) = D 


N N 


J2 Sd) T R-'E{vd)v(j) T )R-'S(j ) 
i=lj=l 


\D (18) 


When the output residuals are assumed to be zero mean white [cf. 
Eq. (5)], then 


(ID 


E[vdMj) T )=RS ij 


(19) 
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From Eqs. (14), (18), and (19), it is easy to see that the parame- 
ter covariance matrix reduces to the dispersion matrix D when the 
output residuals are white. 

By definition, when v is a zero mean weakly stationary random 
process, 13 


with measurement equations 

a m (i) =ct (i) + uj(0 


q m (i) = q(i) + U 2 (0 


E {v(i)v(J) T } = 9t„(* - j) = X„(j - i) (20) 


where — j) is the autocorrelation matrix for the output residual 
vector. The estimated parameter covariance matrix can be computed 
by substituting for E{v(i)v(j) T } from Eq. (20) into Eq. (18) and 
using an estimated value for ~j).An estimate for (i -;) can 

be obtained using the colored residuals from conventional maximum 
likelihood estimation. The autocorrelation estimate accounts for the 
frequency content of the colored residuals in the expression for the 
parameter covariance [Eq. (18)]. Substituting Eq. (20) into Eq. (18) 
results in 


cov (0) 


N N 

= D J>(0 r *“ 1 £ 

_ i ~ 1 j ~ 1 


3 Mi - j)R-'S(j) 


D (21) 


where 9?„(i - j) is computed using the discrete unbiased estimate 
of the output residual autocorrelation 13 : 


1 N ~ k 

&«(*) - JjZTk H v(,)v( ' + k)T = ( 22 ) 


Equation (21) is the expression for the parameter covariance ma- 
trix for colored residuals, which are correlated in time. Equation (22) 
was used to estimate 9Mi - j) in Eq. (21). The values forD, R l , 
and S are from the conventional maximum likelihood estimation. 
Equations (21) and (22) embody the postprocessing applied to a 
conventional maximum likelihood solution to account for colored 
residuals. 

Because this technique postprocesses the output residuals from 
conventional maximum likelihood estimation to correct the Cram6r- 
Rao bounds, all of the properties of conventional maximum likeli- 
hood parameter estimation in a practical flight test data analysis 
application remain unchanged. These properties are discussed in 
Refs. 2 and 3. 

For equation error parameter estimation, the model has a single 
output 

z(i) = x(i) r 0 + e(i) i = l,2,...,// (23) 

where x(i) is an n p x 1 vector of regressors at the ith data point 
and e(i) is the equation error. The preceding analysis applies to this 
case as well, and the equivalents of Eqs. (21) and (22) are 


cov (6 ) = D 


V N 

£*(«) m “ (i ~ 2)*0') T p (24) 

_ 1 = 1 j = 1 J 


and 


where 


and 


= jfZk J 2 e(i)e(i + k) = &«(-*) 


D = 


L-=i 


e(«) = z(i) -x(i) T 0 


(25) 


(26) 

(27) 


Results 

The longitudinal short period dynamics of the F-l 8 HARV fighter 
aircraft at approximately 20-deg angle of attack were studied. The 
model state equations in wind axes are given by 


a= £[ Cz - a+c *<w +c «> s ' +c ‘* 


+ q 


+ |r(cos(<I> m ) cos(© m ) cos(a) + sin(0 m ) sin(a)] (28) 

q = (qS£/l y )[c Ma <x + C Mq (qc/ 2V) + C Mi &, + CJ,J 




o-Hf, 

mg l 


Cz.o(i) + Cz, + Czi'SsO) + a* 


+ vs (0 


1 = 1.2 N (29) 

assuming that Cz ^ and a z & a tw . Initial conditions for the 
states were computed from the measured time histories of a and q 
using a time domain smoother. 14 The parameters Q o , C* M , and a* 
include both aerodynamic and measurement biases. ° 

To validate the new technique for computing Cramer-Rao 
bounds, 200 Monte Carlo simulation runs were made using various 
colored measurement noise processes. Each noise sequence had part 
of its power in the frequencies between 0 and 1 Hz inclusive (roughly 
the frequency band of the uncorrupted simulation outputs), with the 
remaining power taken up by white Gaussian noise out to the Nyquist 
frequency. The narrow-band portion of the colored noise sequence 
was generated by passing zero mean white Gaussian noise through 
a fifth-order Chebyshev low-pass filter with frequency cutoff set at 1 
Hz. The resulting narrow-band noise was combined with wideband 
noise from a separate realization of the zero mean white Gaussian 
noise process. The percentage of the total noise power from the 
narrow-band noise was determined by a random number with uni- 
form distribution on the interval [0, 1 00]. The resulting colored noise 
sequence was then scaled to achieve approximately a 5/1 signal to 
noise ratio for the simulated noisy output This procedure was car- 
ried out for each individual simulated output on each Monte Carlo 
run. Figure 1 shows the power spectral density for the colored noise 
added to cl for run 200, where 19% of the noise power was in the 
frequency range of 0-1 Hz, inclusive. Colored noise sequences gen- 
erated in this way are representative of residual sequences observed 
when analyzing real flight test data and were used for that reason. 

To make the Monte Carlo simulation runs realistic, the stabilator 
input was taken from measured data for the F-l 8 HARV flying a 
maneuver designed specifically for accurate parameter estimation. 15 
The stabilator input is shown as the solid line in Fig. 2. The param- 
eter values used in the simulations (given in column 2 of Table 1) 


Table 1 Parameter estimation results for typical Monte Carlo runs 



0 


Run 47 



Run 185 


e 

n 

r)c 

e 

T) 

r}c 

c z* 

-2.00 

-1.30 

4.24 

0.94 

-3.78 

9.33 

2.26 

c z< 

-65.0 

-13.8 

4.10 

0.99 

-81.8 

1.37 

0.33 


-0.90 

-1.48 

4.20 

0.95 

0.04 

6.46 

1.59 

c l 

-0.80 

-1.08 

4.73 

1.04 

0.01 

10.46 

2.33 

Cm* 

-0.30 

-0.35 

7.00 

1.67 

-0.32 

3.09 

0.48 

CM q 

-16.0 

-17.3 

3.55 

0.86 

-13.2 

7.04 

1.13 

C M S J 

-0.70 

-0.67 

5.42 

1.03 

-0.72 

3.86 

0.57 

C Mo 

0.08 

0.09 

5.77 

1.20 

0.09 

5.01 

0.73 

< 

-0.70 

-1.07 

6.11 

1.36 

0.04 

9.48 

2.14 


0.07 

p 0.06 

v i v i 

0.05 

0.04 

0.03 

0.02 

0.01 

0 


Fig. 1 Example simulated colored noise power spectrum, run 200. 



frequency (hz) 
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Fig. 2 SUbilator input time histories from five repeated flight test ma- 
neuvers. 



Range 

Fig. 3 Cjn a estimates from Monte Carlo simulation. 

approximately reflect the short period dynamics of the F- 18 HARV 
at 20-deg angle of attack. The stabiiator input and parameter values 
were the same for each simulated data am, so that the information 
in the data was constant from run to run. The sampling rate was 
50 Hz, and the data record length was 14 s. Maximum likelihood 
estimation as described in the previous section was used to estimate 
the parameters. 

Since the true parameter values were known for the simulated 
data, the true accuracy of the maximum likelihood estimates could 
be compared to the accuracy indicated by the Cram6r-Rao bound 
calculations. The conventional Cram6r-Rao bounds for the param- 
eter standard errors were denoted by a and were computed from 
Eqs. (14) and (15). The Cram6r-Rao bounds for the parameter stan- 
dard errors corrected for colored residuals were denoted by cr c and 
were computed as the square root of the diagonal elements of the co- 
variance matrix from Eq. (21), using Eq. (22) to estimate the output 
residual autocorrelation. Results from both the conventional com- 
putation and the corrected calculation were expressed in terms of 
the ratio of the absolute deviation of each parameter estimate from 
its true value to the computed Cram6r-Rao bound for the parameter 
standard error. This accuracy measure was assigned the symbol 77 : 

i? m\9 -0 \ M rj c =\e — 0\/a c (30) 

For a maximum likelihood estimator, the probability distribution 
of the parameter estimates approaches a Gaussian distribution cen- 
tered on the true value as the number of data points gets large. 
Evidence of this can be found in Fig. 3, which is a histogram of 
the parameter estimates from all 200 Monte Carlo runs for the Cu„ 
parameter. Corresponding histograms for the other estimated pa- 
rameters were similar. When a equals the standard deviation of the 
population of parameter estimates, the quantity (0 — 0)/a is a stan- 
dardized normal deviate. It follows that 77 < 3 nearly all of the time 
when a is representative of the scatter in the parameter estimates, 
because a standardized normal deviate lies within ±3 standard devi- 
ations of the mean 99.7% of the time. Analogous statements apply 
for 77c and cr c . 


Table 1 shows results for two representative Monte Carlo runs. 
Columns 4 and 5 for run 47 and columns 7 and 8 for nin 1 85 show 
that the corrected Cram 6 r-Rao bounds accurately reflected the true 
parameter accuracy, whereas the conventional Cram 6 r-Rao bounds 
were optimistic (i.e., too small) and produced 77 ratios that exceeded 
3 for almost every estimated parameter. Considering the full set of 
200 Monte Carlo runs, Tfcble 2 gives the mean values and standard 
errors of 77 and t\ c for each estimated parameter. These data show 
that the conventional Cram 6 r-Rao bounds were inaccurate on the 
average and exhibited a large variability, whereas the converse was 
true for the corrected Cram£r-Rao bounds. 

Table 3 gives another summary of the parameter estimation results 
for the 200 Monte Carlo simulation runs. The second column of the 
table gives the mean values of the parameter estimates, and the third 
column gives the sample standard errors for the parameter estimates, 
computed from the scatter of the parameter estimates and denoted 
by s. Columns 4 and 6 give the mean values of the Cram 6 r-Rao 
bounds for the parameter standard errors computed using the con- 
ventional and corrected techniques, 0 and & ct respectively. Columns 
5 and 7 show the ratio of the sample standard errors for the param- 
eter estimates to d and a c , respectively. These values are far less 
than 3 for the corrected calculation of the Cramdr-Rao bounds, in- 
dicating a proper accounting for the changes in the residual spectra, 
whereas the conventional calculation of the Cram 6 r-Rao bound was 
optimistic, producing values of the s/<x ratio greater than 3. 

The data in Tables 1-3 demonstrate that the extent to which 
the conventional Cram 6 r-Rao bounds misrepresented the true 


Table 2 Parameter accuracy statistics 
tor 200 Monte Carlo runs 



Conventional 

Corrected 


<?n 

Vc 

°V 

C Za 

3.84 

2.85 

1.37 

1.10 

Cz q 

3.43 

2.72 

1.10 

0.83 

Czi, 

3.30 

2.51 

1.10 

0.91 

ci 

3.88 

2.89 

1.36 

1.05 

Cm. 

4.03 

3.42 

1.44 

1.38 

Cm, 

3.62 

2.66 

1.26 

1.08 

Cu h 

3.63 

2.72 

1.23 

0.88 

C k 

3.97 

3.19 

1.37 

1.15 

a l 

3.88 

2.80 

1.39 

1.06 


Table 3 Parameter estimation results for 200 Monte Carlo runs 
Simulation Conventional Corrected 

0 s a s/a o c sfa c 

Cza —1.94 0.819 0.179 4.56 0.521 1.57 



frequency (hz) 

Fig. 4 Power spectrum of residuals, flight test run X. 
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parameter accuracy was neither consistent nor predictable from pa- 
rameter to parameter or from run to run. This phenomenon has been 
observed previously when analyzing flight test data from repeated 
maneuvers. 7 It follows that the common practice of applying a fixed 
correction factor to the conventional calculation of the Cram£r-Rao 
bounds is incorrect to a varying and unpredictable degree in cases 
where coloring of the residual spectra varies, as in this simulation 
study. Changes in the coloring of the residuals similar to those stud- 
ied here can easily be brought about in practice by changes in the 
model structure, the maneuver, the flight condition, or the instru- 
mentation. 

Next, flight test data were analyzed from five repeats of the same 
longitudinal maneuver, flown on the F-18 HARV at approximately 
20-deg angle of attack and 25,000 ft. The input was applied to the 



frequency (taz) 

Fig. 5 Power spectrum of q residuals, flight test run 1. 



symmetric stabilator by a computerized onboard excitation system 
(OBES), so that the runs were veiy nearly repeats of one another. 
Figure 2 shows the excellent repeatability using the OBES system 
for five repeated runs of the stabilator input maneuver. All of the 
data used for analysis were sampled at 50 Hz. Corrections were 
applied to the angle-of-attack and accelerometer measurements to 
account for sensor offsets from the center of gravity, and the angle- 
of-attack measurement was corrected for upwash. Data compatibil- 
ity analysis 16 revealed that the data from the sensors were consistent 
to a degree that warranted no further corrections. The same model 
given in Eqs. (28) and (29) was used for the flight test data analysis, 
with measured time histories used for V and q . Maximum likelihood 
parameter estimation was carried out using the procedure described 
in the preceding section. 

Table 4 gives flight test results in a format similar to Tbblt 
3. Column 7 shows that the corrected Cram^r-Rao bounds were 
an accurate measure of the scatter in the parameter estimates. In 
column 5, the conventional Cramer— Rao bounds were again opti- 
mistic for the pitching moment (Cm) parameters but were close to 
correct for the vertical force (C 2 ) parameters. The reason is that the 
a and a z measurements are the main influences on the Cz parame- 
ters, and the residuals for both these outputs exhibited considerable 
power at high frequencies, due to unmodeled structural modes. The 
power spectrum for a typical a z residual (from run 1) is shown in 
Fig. 4. These colored residual spectra roughly resembled constant 
power out to the Nyquist frequency, which is the assumption made 
in the theory underlying the conventional Cram6r-Rao bound cal- 
culation. The q measurement did not have these high-frequency 
components, and so the conventional Cram6r-Rao bound calcula- 
tion gave very optimistic values for the C M parameters. The power 




Fig. 6 Parameter estimates with conventional and corrected error bounds: left bar = conventional error bound and right bar = corrected error 
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Table 4 Parameter estimation results from flight test data 



Right 

Conventional 

Corrected 

e 

s 

a 

s/a 

<r c 

s/a c 

CZa 

-2.04 

0.161 

0.083 

1.95 

0.192 

0.84 

Cz, 

-61.0 

4.86 

3.48 

1.40 

7.67 

0.63 

c Zl , 

-0.92 

0.081 

0.044 

1.82 

0.081 

1.00 

C l 

-0.72 

0.068 

0.031 

2.21 

0.068 

1.00 

Cm, 

-0.30 

0.065 

0.006 

10.80 

0.052 

1.26 

Cm, 

-19.1 

2.55 

0.372 

6.86 

2.84 

0.90 

Cm,, 

-0.74 

0.048 

0.007 

6.90 

0.051 

0.94 


0.08 

0.025 

0.002 

10.62 

0.020 

1.24 

< 

-0.66 

0.067 

0.031 

2.17 

0.070 

0.95 


spectrum for a typical q residual (from run 1) is shown in Fig. 5. The 
corrected calculation of the Cram6r-Rao bound worked equally well 
for the C z and Cm parameters because information about the partic- 
ular coloring of the residuals was incorporated automatically via the 
autocorrelation estimate from Eq. (22) used in Eq. (21). For the C Ma 
parameter, the ratio of parameter estimate scatter to the corrected 
standard error was slightly higher than for the other estimated pa- 
rameters (s/5 c = 1.26 from Table 4, column 7). This anomaly was 
traced to the fact that the control law for the F-l 8 HARV scheduled 
leading- and trai ling-edge flap deflection in proportion to the angle 
of attack. The estimated C Ma and parameters, therefore, in- 
cluded the effects of the flap deflections, which essentially changed 
the wing camber. Run 5 began from a slightly different initial trim 
condition, which caused both leading- and trailing-edge flap deflec- 
tions to differ from the other runs by 1-2 degrees throughout the 
maneuver. The change in the effective wing camber was manifested 
in the C Ma parameter estimate, but not the estimate, in agree- 
ment with slender wing theory and wind-tunnel data. 17 The true Cm„ 
parameter, therefore, contained variation due to the changes in the 
effective wing camber, violating the assumption that the parameters 
should be constant for repeated maneuvers. The result was a larger 
scatter in the estimates of C Ma , which increased the value of s/d c . 
A similar effect was seen for C* Mo , whose estimated value reflects 
the trim condition to some extent 

Figure 6 depicts the parameter estimation results for the aerody- 
namic parameters. The error bars to the left of the round symbols 
marking the individual parameter estimates represent the Cram6r- 
Rao bounds for the standard errors computed using the conventional 
calculation The error bars to the right of the parameter es- 

timate symbols represent the Cram6r-Rao bounds for the standard 
errors from the corrected calculation (±la c ). These plots and the 
accompanying data in Table 4 show that the standard calculation 
for the Cram6r-Rao bounds gave optimistic values compared to the 
scatter in the estimates from repeated maneuvers, whereas the cor- 
rected calculation for the Cramdr-Rao bounds produced Cramdr- 
Rao bounds that accurately reflected the scatter of the estimates. 

Concluding Remarks 

Algorithms for aircraft parameter estimation using the output er- 
ror formulation of maximum likelihood are in widespread use. The 
Cramdr-Rao bounds characterizing parameter accuracy that are ob- 
tained from conventional calculations are known to be generally op- 
timistic (i.e., too small) in practice, compared to the scatter in param- 
eter estimates from repeated maneuvers. Estimated parameters have 
limited utility when there is no firm idea of their accuracy. In this 
work, an expression for correcting Cram6r-Rao bounds from max- 
imum likelihood estimation with colored residuals was developed 
and validated. This result is important because the residuals from 
maximum likelihood estimation are almost always colored in prac- 
tice, due to deterministic modeling error. The technique was shown 
to be applicable to equation error parameter estimation as well. 

The calculations involved in the algorithm for computing 
Cram6r-Rao bounds that account for colored residuals can be 
carried out in a short subroutine called at the conclusion of a 


conventional maximum likelihood estimation algorithm. Bandwidth 
of the dominant power in the residuals need not be known or es- 
timated because it is accounted for automatically in the algorithm 
by an unbiased estimate of the residual autocorrelation. There is no 
need for correction factors. The algorithm was shown to work for 
a wide range of colored residual spectra similar to what might be 
encountered in real flight test data analysis. All calculations are per- 
formed in the time domain, obviating the need for frequency domain 
analysis of the residuals. 

The corrected calculation for the Cram6r-Rao bounds presented 
here produced consistently accurate measures of the scatter in the pa- 
rameter estimates, using an algorithm with moderate computational 
cost that was applied as a postprocessing of the output residuals 
from a conventional maximum likelihood solution. 

Monte Carlo simulation runs using various colored noise se- 
quences were carried out to validate the algorithm. Analysis of flight 
data from repeated maneuvers flown on the F-l 8 HARV demon- 
strated the validity of the technique for computing appropriate pa- 
rameter accuracy measures using real flight test data. The algorithm 
described in this work was shown to be an effective means of ac- 
curately determining the quality of parameter estimates from the 
output error formulation of maximum likelihood estimation. 
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